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INTRODUCTION 


Data insufficiency, poorly conditioned matrices and singularities in 
equations occur regularly in complex optimization, correlation, and 
interdisciplinary model studies. This work concerns itself with two methods of 
obtaining certain physically realistic solutions to ill-conditioned or singular 
algebraic systems of linear equations arising from such studies. 

Two efficient computational solution procedures that generally lead to 
locally unique solutions are presented when there is insufficient data to 
completely define the model, or a least-squares error formulation of this system 
results in an ill-conditioned system of equations. 

If it is assumed that a "reasonable" estimate of the uncertain data is 
available in both cases cited above, then we shall show how to obtain realistic 
solutions efficiently, in spite of the insufficiency of independent data. 

The proposed methods of solution are more efficient than singular-value 
decomposition [1] for dealing with such systems, since they do not require 
solutions for all the nonzero eigenvalues of the coefficient matrix. 


OBJECTIVES 


REVIEW MATHEMATICAL FORM OF ILL-CONDITIONING 

OUTLINE PHYSICAL SOURCES OF ILL-CONDITIONING 

REVIEW SINGULAR-VALUE DECOMPOSITION SOLUTION APPROACH 

PRESENT A COMPUTATIONALLY SIMPLER METHOD 

OFFER MATHEMATICAL & GEOMETRIC INTERPRETATION OF SOLUTIONS 

DISCUSS PHYSICAL INTERPRETATION AND RELATIONSHIP TO: 

1. MATH MODEL IMPROVEMENT 

2. GRADIENT METHODS 

3. BAYESIAN ESTIMATION 
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PROBLEM STATEMENT 


Given the physical system with "m" bits of experimental data, j Y g j and "n" 
uncertain parameters $ r } , * 

assume that a math model exists which yields the "m" values f Y I based 
upon the "m" estimated parameters | r ^ . *■ 

, s It is desired to determine the "best" fit of parameters [r} so that 
Jy - y|J i s a minimum. 

A one-term Taylor approximation yields: 



where is the Taylor series truncation error plus experimental error. 
The mXn sensitivity matrix [S] is then defined below. 



, J { y a | 


IS1 * 1 A| 

SENSITIVITY MATRIX: 

M X N [HdRi _ 

NO. PARAMETERS: 

{A R i = 1 R - R 0 1 


N X 1 

SOLUTION AND OUTPUT DATA: 

{AY} = f Y e - Y a } 


M X 1 

RESIDUAL: 

{ R { = {AY( - IS1 [a R } 


M X 1 

TAYLOR SERIES: 

[Y E j = { Y a ] ♦ I S I f 4 r f ♦ Res 


LEAST-SQUARES MINIMUM RESIDUAL: IS T S1 U r } = IS1 T fd Y} 
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CONDITIONING 


1. If m>n, there are more data than parameters, and so we shall seek a least- 
error-squared fit. 

2. If m< n, then there are an infinite number of solutions and it is not clear 
which one to use. 

T 

Let us pursue case (1) first (i.e. m >n) . If [S S] is numerically 
invertible, then a least-error-squared solution is possible by efficient 
triangular factorization of [S^S] into [LL^] , and subsequent forward and 
backward substitution. However, if "n" is only moderately large, then solutions 
can become difficult due to ill-conditioning. Also, if the rank of [S S] is 
less than n, factorization is not possible. Now, if m<n and a solution is 
somehow obtained, what is its interpretation? 

One possible approach to these problems is to use the singular-value 
decomposition (SVD) solution approach [1]. This technique requires obtaining all 
the q (in) nonzero eigenvalues ['^s] of [S T S] and the corresponding nxq modal 
matrix, [U] . If the rank of [S T S] is n, then the interpretation of SVD is 
clear. However, if the rank is less than n (as it may be for m>n and will be 
for m^n), interpretations of the SVD solutions are not obvious. 


M > N 

: 1 . 

S T S 

WELL CONDITIONED 

(NUMERICALLY & PHYSICALLY) 


2. 

S T S 

ILL-CONDITIONED 

(NUMERICALLY, BUT NOT PHYSICALLY) 

M < N 

: 3. 

S T S 

ILL-CONDITIONED 

(NUMERICALLY & PHYSICALLY) 


CASE OF WELL-CONDITIONED NUMERICALLY, BUT NOT PHYSICALLY, IS NOT CONSIDERED 


IF 

CASE 1: 

IS T S] 

= (LI 

IU T , L = CHOLESKY DECOMPOSITION = IV;*. 0 I 

IF 

CASE 2: 

is T sf‘ 

Jw U 

I'A.]" 1 U T SINGULAR-VALUE DECOMP. 



N X N 

NXQ 

QXQ QXN 



X = (Q 

£ N) 

NON-ZERO EIGENVALUES OF IS T SI 


IUKORRESPONDING MODAL MATRIX 
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EPSILON DECOMPOSITION 


The previously discussed procedure is laborious at best. A much simpler 
procedure is now presented. First, pick a small value £ , and set 

[ S T S ] = [S T S] + <£ [I] 

T T 

such that [S S 1 may be efficiently factored into [LL ] where [L] is lower 

triangular. This is then used to solve 
[S3] (Ar\ = [S] T f Ay} 

Next, decrease £ until {ir| approaches an asymptote, or the factorization of 
[S S] into [LLj] breaks down due to ill-conditioning. 

This method has been termed a "Levenberg-Marquardt" method by Luenberger 
[2], who considers it to be a modification of Newton's Method and steepest- 
descent (see pp. 225-227 of Reference 2). 


( SIMPLER PROCEDURE ) 


IS 


T 


SI 


is T si ♦ £ m = i l i iu T 

£>£>£>•><’ 0 

*12 3 


BUT, WHAT DO SOLUTIONS MEAN? 


EXPLANATION BASED ON SIMPLE 2-DOF INTERPRETATION 
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TWO DOF EXAMPLE 


or 


Consider the two degree-of-freedom (DOF) example with m=l, where we are given 
r S 3 = [a b] 

Therefore, 


T 

[S A S] = 


T - 1 
[S T S] = 


a 2 + £ 


a*b 


a *b 

b 2 +£ 


S T S | = 


£ (£ + a 2 + b 2 ) 


b 2 +£ 


-a b 


0 


-a b 


a 4 +£ 


. ij T f 7 

Thus, the least-error-squared solution becomesjArj = [S S] [S] )Ay ( 

a Ay 


lim (Arj 

g—O 


lim 

e— O 


a 2 + b 2 + £ 
b Ay 
a 2 + b 2 +£ 


a Ay 

2 . ,2 

a + b 

b 4 y 

2 . ,2 
a + b 
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GEOMETRIC INTERPRETATION 


A geometrical interpretation of the solution shows that the ideal solution 
occurs for the point on the line [S]{Ar| y \ which is closest to the origin. 

Note that if £ > 0, then the solutions fall short of this ideal point. 

A comparison of our proposed epsilon-decomposition approach with the method 
of singular-value decomposition yields the identical result for[dr j , but 
through different computational steps. 



1245 



mathematical interpretition 

The previous two-pages suggest that an n-DOF solution approach is not 
only to minimize the sum of the residuals squared but also to minimize the distance 
on this line to the origin of{Ar} space, i.e. 

x T T 

<J> = R R +<fdr 4r 

T T 

If [S S] is not ill-conditioned, simply set £ to zero and solve by [LL ] 

decomposition. However, if it is ill-conditioned, choose an £ and solve by [LL ] 
decomposition. Keep reducing £ until the solution asymptotes to certain values, 
or this type of factorization breaks down numerically. 

Choosing a starting <£ is arbitrary, but the two-DOF ^xample suggests that £ 
should be small compared to the smallest eigenvalue of [S S]. However, it is 
dangerous to generalize. 

We have shown that this procedure works, not only when m>n, but also when 
m <n (which was initially called case 2 here). However, a more direct solution 
is possible if m<n and this is shown on the following page. 


f R ] ={4 Y ] - (SI {A R j 

4> = { R } [R J +£ [4 Rj T [4 R | 


- f 0 } x ((S T SJ ♦ £ m>f4R J = IS] J [A Y / 




z i 

^)fdR} 
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UNDER-DETERMINED SYSTEMS 


—For the case where m<n, [S S] will not exist. However, if the rank of 
[ S rj.S ] is m, then its inverse will exist, and it is possible to show that 
[S S] -1 [S] T is equivalent to S T [SS^]' (as £ approaches zero). Therefore, it is 
not necessary to use epsilon in determining { A r } , for m<n, when the rank of [SS T J 
is m. However, if its rank is less than m, we may still wish to use the 
present method on either [S^S] or [ SS x ] . Our method will yield almost identical 
results if epsilon is sufficiently, small, as compared to the smallest non-zero 
eigenvalue of either [S^S] or [SS T ]. 


M X M 

M<N .Case.. l i -S ll JS £L 


Given: [S t S 1 f a r } ■ [S1 t {4Y) lS T si * 0, [ar 1 * 

/ NXN 

Consider: lim IS T S1 [S T 1 = lim (S t S+£I) -1 [S T ] 

£—-0 &~~o 

[S T S] [S) T = LIM <S T S *£ I) S T 

*— +-0 

Assume: I S S T I=* 0, Post-multipiy Eq (1) by IS S T )' 1 



.'.S T 

= LIM (S T S 

e—-o 

+ £ tin s T is s 1 r 1 

Assume: 

l S T S 

1*0, e. > 0 
Premult. Eq(2) 

BY LIM [S T S ] _1 

s-—0 

Then. 

LIM 

£ 

[$ T S ) _1 S T 

0 

* S T IS S T l " 1 

and 

SO 

[AR\ 

S T IS S T J " 1 [a y} 


tS T Sl _1 tSJ T fAY} 


( 1 ) 


, . ( 2 ) 
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PHYSICAL INTERPRETATION OF PROPOSED PROCEDURE 


Strictly speaking, the unique solution of an ill-conditioned system of 
equations is nonsense. What we have shown is that a unique solution is possible 
if an additional condition is introduced, in the form of a constraint, upon the 
minimization problem; namely, that the solution must not only satisfy the 
original "m" equations in a least-squared error sense, but if the resulting nxn 
system is ill-conditioned, then the proper solution should be closest to the 
initial math model. This latter condition is sufficient to uniquely determine 
the solution in most cases. In our minimization formulation, epsilon represents 
a weighting of the closeness-to-the-origin constraint, Ir-r^ | = min, 
relative to the size the m-equation residual, IRJ , 

minimization. Thus, if we wish to minimize I Rl € should be made as small as 
practically possible while still delivering a unique solution forfdr } . 



(NON-UNIQUE SOLUTIONS TO S , 8 S a EXIST ALONG INTERSECTION ) 

$ “ { R 1 ♦ £ fW R) ■ Min (MAKES PROBLEM WELL-POSED 8 SOLUTION UNIQUE) 

[*R I - { 0 / IS THE INITIAL GUESS 
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RELATIONSHIP TO EXISTING METHODS 


Physically speaking. our approach implies that the initial solution-guess, 

J r 0 ^, is reasonable and that if the formulation leads to some insufficiency in 
uniquely defining the system (i.e. an ill-conditioned [S^S] matrix) then the 
criterion of closeness of { r J to {r Q \ shall be imposed to uniquely define the 
system. By making £ as small as is practical, we are simply weighting this 
closeness of { rj to { r Q } condition as secondary to the least-square-error 
criterion. Viewed in this way, we may consider a different weighting of the 
various parameters of { r } through the diagonal weighting matrix O W r *,]. 
Continuing this concept for the various residuals, [ R| , as well as through the 
weighting matrix [Wy] we arrive at a Bayesian formulation i.e. to find the 
minimum of where * is given by: 

$ = [R> T [W y ] [ R } + {r - r 0 } T ['W ? ] (r - r Q j . 


NEWTON METHODS: 


> f 2 TRAVEL ALONG GRADIENT TO 

SOLUTION FROM INITIAL GUESS 


BAYESIAN ESTIMATION: 
MINIMIZE: 




{ R } T [W y 1 { R { 

t 

Weighted 

Residuals 


{dRl T IW R l f/iR j 

t 

£ (I) 

Proposed Solution 
Approach 
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CONCLUSIONS 


We have derived a numerical technique for solving ill-conditioned systems of 
equations which does not rely upon computation of the system’s non-zero 
eigenvalues and eigenvectors (as is necessary with singular-value 
decomposition) . While others have proposed this same procedure, it is felt that 
we have given it a more formal and less heuristic derivation. In addition, we 
have supplied a physical interpretation which gives insight into the results 
which the method yields, as well as its relation to Bayesian estimation methods. 


AVOIDED USE OF SINGULAR-VALUE DECOMPOSITION 
(COMPUTATIONALLY HARDER) 

PROVIDED MATHEMATICAL BASIS FOR PROPOSED SOLUTION APPROACH 

PRESENTED A GEOMETRIC INTERPRETATION 

PRESENTED A PHYSICAL INTERPRETATION 

IMPROVED INSIGHT INTO SOLUTION OF ILL-CONDITIONED EQS, 
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SYMBOLS 


f r ! 

i*.\ 

{Ar| 

W 

w 

[ay] 


£ 

[I] 

[L] 

[S] 

[U] 

[Vr] 

[Wy] 

[\M 

T 

-1 

-P 


Elements of a matrix 

Parameters (or solution) of math model being sought 
Starting estimate of math model parameters 
r - r 

o 

Analytically determined solution vector 

Experimental determined values 

Y - Y 
e a 

Residual 

2-dimensional surfaces in 3-dimensional space 

Epsilon parameter 

Functions to be minimized 

Identity matrix 

Lower triangular matrix 

Senitivity matrix 

Modal matrix corresponding to non-zero eigenvalues 
Initial parameter confidence weighting matrix 
Experimental data confidence weighting matrix 
Diagonal eigenvalue matrix 

Transpose of a matrix (used as a superscript) 
Inverse of matrix (used as a superscript) 
Pseudo-inverse (used as a superscript) 
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